\documentclass[a4paper,11pt]{article}
\usepackage[top=3.5cm, bottom=4cm]{geometry}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{lmodern}
\usepackage[francais]{babel}
\usepackage{hyperref}
\usepackage{pgf}
\usepackage{tikz}
\usepackage{pgfplots}

\usepackage{pdflscape}
\usetikzlibrary{arrows,automata}
\usepackage{amssymb}
\hypersetup{urlcolor=blue,linkcolor=black,citecolor=black,colorlinks=true} 
\usepackage[french,boxruled,inoutnumbered]{algorithm2e}	
\usepackage{fancyvrb}
\usepackage{graphicx}
\usepackage{float}

\fvset{
frame=single,
fontsize==\footnotesize , 
%numbers=left,
}

%%%% debut macro pour enlever le nom chapitre %%%%
\makeatletter
\def\@makechapterhead#1{
  \vspace*{0\p@}
          {\parindent \z@ \raggedright \normalfont
            \interlinepenalty\@M
            \ifnum \c@secnumdepth >\m@ne
            \Huge\bfseries \thechapter\quad
            \fi
            \Huge \bfseries #1\par\nobreak
            \vskip 40\p@
}}

\def\@makeschapterhead#1{
  \vspace*{50\p@}
          {\parindent \z@ \raggedright
            \normalfont
            \interlinepenalty\@M
            \Huge \bfseries  #1\par\nobreak
            \vskip 40\p@
}}
\makeatother
%%%% fin macro %%%%

\title{CUDA\\Fractale Julia}
\author{Allan BLANCHARD, Sébastien POLNECQ}

\begin{document}

%%PAGE TITRE
\begin{titlepage}
\begin{center}
\includegraphics[scale=0.6]{polytech.png}
\rule[1ex]{\textwidth}{0.1mm}
{\huge
  CUDA\\
  Fractale Julia
}
\rule{\textwidth}{0.1mm}
\vfill
\begin{minipage}{\textwidth}
\centering
{\Large \textsc{Allan BLANCHARD, Sébastien POLNECQ}}
\end{minipage}

\vfill

{\large \today}

\end{center}
\end{titlepage}

\nocite{*}

\section{But Projet Cuda}
\paragraph{}Le but du projet est d'optimiser et d'améliorer les performances système concernant le calcul d'une fractale de Julia.
\section{Julia Init}
\paragraph{}Pour cela, nous sommes partis du code initial fournis qui permettait de lancer un thread par bloc. Chaque bloc représentant un pixel. Les performances du code initial étaient les suivantes : 
\begin{itemize}
\item time kernel = 2057ms
\item time copy = 280ms.
\end{itemize}
\paragraph{}L'ensemble des calculs de performances sera réalisé avec une image de dimension 8192.
\section{Julia Amélioration 1 : pensé Cuda basique}
\paragraph{}En du code de base, nous avons écris une première ``optimisation'' permettant de lancer 256 threads par bloc, les threads étant répartis selon une grille 16x16. 
\paragraph{}Chaque thread de chaque bloc représente également un pixel de l'image. La version précédente ne prenait pas en compte le fait qu'un bloc (physique) possède plusieurs threads physiques ; un grand nombre d'entre eux n'étaient donc pas utilisés.
\paragraph{}Les performances du programme suite à cette prise en compte sont : 
\begin{itemize}
\item time kernel = 35.7ms
\item time copy = 280ms.
\end{itemize}
\paragraph{}Nouvel appel kernel :
\begin{Verbatim}
  #define BLOCK 16
  
  //...

  dim3 grid(size/BLOCK, size/BLOCK);
  dim3 block(BLOCK, BLOCK);
  kernel<<<grid,block>>>(dev_image, size);
\end{Verbatim}
\paragraph{}Dans le kernel on peut voir que le calcul du pixel traité diffère beaucoup. En effet le thread doit se placer à la fois en fonction de la position du bloc dans l'image et en fonction de sa position propre dans le bloc.
\begin{Verbatim}
__global__ void kernel(unsigned char* image, int size) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  int j = blockIdx.y * blockDim.y + threadIdx.y ;

  unsigned int offset = (j*size + i) * 4;

  image[offset] = 0;
  image[offset+1] = 0;
  image[offset+2] = 255*julia(i, j, size);
  image[offset+3] = 255;
}
\end{Verbatim}
\section{Julia Amélioration 2 : Déroulage}
\paragraph{}Cette seconde amélioration consiste à dérouler les boucles de façon à limiter le nombre de sauts et de conditions.
\paragraph{}Nous avons déroulé la boucle concernant le test de l'appartenance d'un complexe au domaine de julia. Nous avons testé avec les déroulages suivants : 2, 4, 8 et 16.
\begin{Verbatim}
for(unsigned int i = 0 ; i < 200/taillederoulage; ++i) {
  //Ce code
  a = a * a + c;
  if(a.magnitude() > 1000) { return 0; }
  //répété 'tailledéroulage' fois
}
\end{Verbatim}
\paragraph{}Les résultats obtenus sont les suivants :
\begin{figure}[H]
  \centering
  \resizebox{8cm}{!}{
    \begin{tikzpicture}
      \begin{axis}[grid=major,
          xlabel=$Optimisation$,
          ylabel=$Temps (ms)$]
        \addplot[smooth,mark=x,red] plot coordinates {
          (1, 35.7)
          (2, 31.0)
          (4, 27.5)
          (8, 25.6)
          (16, 27.6)
        };
      \end{axis}
    \end{tikzpicture}
  }
  \caption{Comparaison déroulage Julia}
\end{figure}
\paragraph{}On constate que le déroulage le plus efficace pour cette fonction est celui de taille 8.
\newpage
\paragraph{}On cherche ensuite à savoir combien de données on doit traiter par threads, de façon à être le plus efficace possible. Le code kernel devient alors :
\begin{Verbatim}
#define STEP 'N'

//...

__global__ void kernel(unsigned char* image, int size) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  int j = (blockIdx.y * blockDim.y + threadIdx.y)*STEP;

  #pragma unroll
  for(unsigned int k(0); k < STEP; k++){
    unsigned int offset = ((j+k)*size + i) * 4;

    image[offset] = 0;
    image[offset+1] = 0;
    image[offset+2] = 255*julia(i, j, size);
    image[offset+3] = 255;
  }
}
\end{Verbatim}
\paragraph{}Les résultats obtenus sont les suivants :
\begin{figure}[H]
  \centering
  \resizebox{8cm}{!}{
    \begin{tikzpicture}
      \begin{axis}[grid=major,
          xlabel=$Optimisation$,
          ylabel=$Temps (ms)$]
        \addplot[smooth,mark=x,red] plot coordinates {
          (1, 25.6)
          (2, 26.8)
          (4, 24.6)
          (8, 29.6)
        };
      \end{axis}
    \end{tikzpicture}
  }
  \caption{Comparaison Déroulage kernel}
\end{figure}
\paragraph{}On remarque que traiter 4 pixels par thread est un bon compromis.

\section{Julia Amélioration 3 : ``Shared''}
\paragraph{}Le but est de réaliser les opérations de calcul dans une mémoire locale et d'envoyer uniquement les résultats au niveau de la mémoire partagée une fois ceux-ci terminés.
\paragraph{}On réserve donc une zone de mémoire locale ({\it shared}) qui représente le bloc que nous allons calculer. Une fois la zone complétement remplie on l'envoie dans la mémoire globale. La réinterprétation des deux tableaux (shared et global) en tant qu'entiers, permet d'améliorer les performances par rapport à un envoi en tant que {\it unsigned char}.
\paragraph{}Le code du kernel devient alors :
\begin{Verbatim}
__global__ void kernel(unsigned char* image, unsigned int size) {
  unsigned int x = blockIdx.x*blockDim.x+threadIdx.x;
  unsigned int y = (blockIdx.y*blockDim.y+threadIdx.y)*STEP;

  __shared__ unsigned char pixels[BLOCK*STEP][BLOCK][4];

  #pragma unroll
  for(unsigned int k = 0; k < STEP; ++k){
    pixels[threadIdx.y*STEP+k][threadIdx.x][0] = 0;
    pixels[threadIdx.y*STEP+k][threadIdx.x][1] = 0;
    pixels[threadIdx.y*STEP+k][threadIdx.x][2] = 255*julia(x,y+k,size);
    pixels[threadIdx.y*STEP+k][threadIdx.x][3] = 255;
  }

  #pragma unroll
  for(unsigned int k = 0; k < STEP; ++k){
    unsigned int offset = ((y+k)*size+x);
    reinterpret_cast<int*>(image)[offset]
      = reinterpret_cast<int*>(pixels[threadIdx.y*STEP+k])[threadIdx.x];

  }
}
\end{Verbatim}
\begin{figure}[H]
  \centering
  \begin{tabular}{| l | c | c |}
    \hline
    Taille image & Sans Shared & Avec Shared\\ 
    \hline
    2048  &  2.2 ms & 1.6 ms\\ 
    \hline
    4096  &  6.5 ms & 6 ms\\ 
    \hline
    8192  &  25.5 ms & 22.8 ms\\ 
    \hline
  \end{tabular}
  \caption{Comparaison avec et sans Shared}
\end{figure}

\section{Julia Bilan Intermédiaire : Amélioration Kernel}
\paragraph{}Graphique intermédiaire des performances en fonction des améliorations apportées :\\
\newline
\begin{figure}
  \centering
  \resizebox{8cm}{!}{
    \begin{tikzpicture}
      \begin{axis}[grid=major,
          xlabel=$Optimisation$,
          ylabel=$Temps (ms)$]
        \addplot[smooth,mark=x,red] plot coordinates {
          (1, 35.7)
          (2, 25.5)
          (3, 22.8)
        };
      \end{axis}
    \end{tikzpicture}
  }
\caption{Comparaison des optimisations}
\end{figure}
\begin{enumerate}
\item Base
\item Déroulage
\item Shared
\end{enumerate}
\paragraph{}Si l'évolution du temps calcul est acceptable, le poids du temps de chargement est, lui, assez élevé par rapport au temps total.
\begin{center}
  \includegraphics[scale=0.5]{tableauinter.png}
\end{center}
\paragraph{}Il nous paraît difficile à partir d'ici d'améliorer le temps kernel. En revanche, il est possible en parallélisant les copies de recouvrir le temps de calcul.
\paragraph{}Dans les parties suivantes, nous nous intéresserons à réduire l'impact de la copie en la réalisant au maximum en même temps que le calcul proprement dit.
\section{Julia Amélioration 4 : ``Streams''}
\paragraph{}Il est possible d'utiliser de la mémoire appelée ``page-locked memory'' avec les streams ; cette zone de mémoire a pour particularité de ne pas être gérée en tant que mémoire virtuelle par l'OS, on a donc l'adresse physique dans la mémoire. Les transferts pourront être directement effectués de la mémoire de la carte graphique vers la mémoire de l'hôte.
\subsection{Sans page-locked memory}
\paragraph{}Le but ici est d'envoyer les ``morceaux'' traités au fur et à mesure de leur création. Par l'utilisation des {\it streams}, il est possible de continuer à effectuer des calculs sur la carte pendant qu'on envoie des données vers le CPU (ou inversement).
\paragraph{}D'après nos tableaux, nous avons considéré que le calcul représentait environ un huitième du temps de transfert mémoire. Nous avons décidé de créer 8 {\it streams} pour améliorer notre programme.
\paragraph{}Les codes suivants ont été modifiés au niveau du code principal :
\begin{Verbatim}
  #define STREAMS 8

  //...

  unsigned int size = static_cast<unsigned int>(atoi(argv[1]));
  size = ((size/(BLOCK*STEP*STREAMS))+(size%(BLOCK*STEP*STREAMS) != 0))
    *BLOCK*STEP*STREAMS;
  
  // On calcule la hauteur d'un bloc donnée de streams
  unsigned int hblock = size/STREAMS;

  //maintenant la hauteur de la grille dépend du bloc stream
  dim3 grid(size/BLOCK, hblock/(BLOCK*STEP));
  dim3 block(BLOCK, BLOCK);

  //....

  unsigned char* dev_image[STREAMS];
  cudaStream_t streams[STREAMS];

  for(uint i(0); i < STREAMS; ++i){
    cudaStreamCreate(&streams[i]);
    cudaMalloc(&dev_image[i], hblock*size*4);
  }
 
  for(uint i(0); i < STREAMS; ++i)
    //on passe la hauteur en cours
    kernel<<<grid,block,0,streams[i]>>>(dev_image[i], size, i*hblock);
  
  for(uint i(0); i < STREAMS; ++i)
    cudaMemcpyAsync(&image[i*hblock*size*4], 
		    dev_image[i], hblock*size*4, 
		    cudaMemcpyDeviceToHost, streams[i]);

  //...
  for(uint i(0); i < STREAMS; ++i){
    cudaStreamDestroy(streams[i]);
    cudaFree(dev_image[i]);
  }  
\end{Verbatim}
\paragraph{}Le kernel a évolué de cette manière :
\begin{Verbatim}
//h représente la 'hauteur' à laquelle se trouve le bloc dans l'image
__global__ void kernel(unsigned char* image, unsigned int size, int h) {
  //pas de changement ici : on agit sur le bloc de donnée 
  //réservé au stream en cours
  unsigned int x = blockIdx.x*blockDim.x+threadIdx.x;
  unsigned int y = (blockIdx.y*blockDim.y+threadIdx.y)*STEP;

  __shared__ unsigned char pixels[BLOCK*STEP][BLOCK][4];

#pragma unroll
  for(unsigned int k = 0; k < STEP; ++k){
    pixels[threadIdx.y*STEP+k][threadIdx.x][0] = 0;
    pixels[threadIdx.y*STEP+k][threadIdx.x][1] = 0;
    //ici en revanche il faut ajouter le déplacement vertical
    pixels[threadIdx.y*STEP+k][threadIdx.x][2] = 255*julia(x,y+k+h,size);
    pixels[threadIdx.y*STEP+k][threadIdx.x][3] = 255;//*/
  }
  // ...
  // pas de changement dans la suite
}
\end{Verbatim}
\paragraph{}Voici les nouveaux temps (totaux), comparés aux anciens :
\begin{figure}[H]
  \centering
  \begin{tabular}{| l | c | c |}
    \hline
    Taille image & Sans Stream & Avec Stream\\ 
    \hline
    2048  &  299.8 ms & 277 ms\\ 
    \hline
    4096  &  76.5 ms & 72 ms\\ 
    \hline
    8192  &  19.6 ms & 21 ms\\ 
    \hline
  \end{tabular}
  \caption{Comparaison avec et sans streams}
\end{figure}
\subsection{Avec page-locked memory}
\paragraph{}On change l'allocation de la mémoire par le biais d'une instruction spécifique CUDA : {\tt cudaMallocHost}. Les {\tt new} d'allocation du tableau sur la machine sont remplacés par le {\tt cudaMallocHost} et les {\tt delete[]} par {\tt cudaFreeHost}.
\paragraph{}L'évolution est la suivante :
\begin{figure}[H]
  \centering
  \begin{tabular}{| l | c | c | c |}
    \hline
    Taille image & No Stream & Stream (no pl) & Stream (pl)\\ 
    \hline
    2048  &  299.8 ms & 277 ms & 120 ms\\ 
    \hline
    4096  &  76.5 ms & 72 ms & 30 ms\\ 
    \hline
    8192  &  19.6 ms & 21 ms & 8 ms\\ 
    \hline
  \end{tabular}
  \caption{Comparaison globale}
\end{figure}
\section{Julia Amélioration 5 : Multicartes}
\paragraph{}Nous avons pour cette amélioration "découpé" le traitement de la fractale en deux: sur les deux cartes à notre disposition. On execute la moitié des Streams sur chacune des cartes.
\paragraph{}Cela permet de réduire le temps de traitement d'un facteur deux.
\begin{Verbatim}
for(uint i(0); i < STREAMS; i++){
  cudaSetDevice(i%2);
  cudaMalloc(&dev_image[i], hblock*size*4);
  cudaStreamCreate(&streams[i]);
}

cudaSetDevice(0);
cudaEventRecord(start0, 0);

cudaSetDevice(1);
cudaEventRecord(start1, 0);

for(uint i(0); i < STREAMS; ++i){
  cudaSetDevice(i%2);
  kernel<<<grid,block,0,streams[i]>>>(dev_image[i], size, i*hblock);
}

for(uint i(0); i < STREAMS; ++i){
  cudaSetDevice(i%2);
  cudaMemcpyAsync(&image[i*hblock*size*4], 
                  dev_image[i], hblock*size*4, 
	          cudaMemcpyDeviceToHost, streams[i]);
}

//...

for(uint i(0); i < STREAMS; ++i){
  cudaSetDevice(i%2);
  cudaStreamDestroy(streams[i]);
  cudaFree(dev_image[i]);
}

//...

cudaFreeHost(image);
\end{Verbatim}
\paragraph{}Ce qui donne les performances suivantes :
\begin{figure}[H]
  \centering
  \begin{tabular}{| l | c | c | c | c |}
    \hline
    Taille & No Stream & Stream (no pl) & Stream (pl) & Multi-GPU\\ 
    \hline
    8192  &  299.8 ms & 277 ms & 120 ms & 61 ms\\ 
    \hline
    4096  &  76.5 ms & 72 ms & 30 ms & 15 ms\\ 
    \hline
    2048  &  19.6 ms & 21 ms & 8 ms & 4 ms\\ 
    \hline
  \end{tabular}
  \caption{Comparaison globale}
\end{figure}
\section{Julia Bilan : "Tableau comparatif"}
\paragraph{}Voici le graphique final des performances (en temps total) pour une image de 8192 (on ne considère pas la version de base parce qu'elle rendrait le graphique illisible, le temps d'exécution étant de plus de 2300 ms) :
\begin{center}
  \resizebox{10cm}{!}{
    \begin{tikzpicture}
      \begin{axis}[grid=major,
          xlabel=$Optimisation$,
          ylabel=$Temps (ms)$]
        \addplot[mark=x,red] plot coordinates {
          (1, 315.7)
          (2, 300.5)
          (3, 300.5)
          (4, 120)
          (5, 61)
        };
      \end{axis}
    \end{tikzpicture}
  }
\end{center}
\begin{enumerate}
\item Base
\item Déroulages
\item Shared
\item Streams
\item Bi-GPU
\end{enumerate}
\end{document}

